#include <iostream>
#include <stdlib.h>
#include <algorithm>

#include "libmri.h"

using namespace std;

int main (int argc, char **argv)
{
	bool adjoint = 0;

	if (argc > 1){
		if (!strcmp (argv[1], "adjoint")) adjoint = true;
	}

	cout << "adjoint? " << (adjoint ? "yes" : "no") << endl;

	shm_array_t *shm_input = 0, *shm_kernel = 0;
	shm_array_t size_output;
	char *output = 0;

	if (0 == shm_input){
		if (get_shm_array("input", &shm_input)){
			perror("input");
			return -1;
		}
	}

	if (0 == shm_kernel){
		if (get_shm_array("kernel", &shm_kernel)){
			perror("kernel");
			return -1;
		}
	}

	if (0 == output){
		output = strdup("output");
	}

	int xs1 = shm_input->dims[0], xs2 = shm_input->dims[1], 
	    xch = shm_input->dims[2], xbatch = 1;
	for (int i = 3; i < shm_input->ndims; ++i){
		xbatch *= shm_input->dims[i];
	}

	int ks1 = shm_kernel->dims[0], ks2 = shm_kernel->dims[1],
	    kch1 = shm_kernel->dims[2], kch2 = shm_kernel->dims[3],
	    kbatch = 1;
	for (int i = 4; i < shm_kernel->ndims; ++i){
		kbatch *= shm_kernel->dims[i];
	}

	if (xbatch != kbatch){
		cerr << "X batch ~= K batch!" << endl;
		return -1;
	}
	int batch = xbatch;

	if ((adjoint and (xch != kch2)) or (not adjoint && (xch != kch1))){
		cerr << "(X,K,adj) Size mismatch!" << endl;	
		return -1;
	}

	memset (&size_output, 0, sizeof (size_output));
	size_output.ndims = shm_input->ndims;
	int ych;
	for (int i = 0; i < shm_input->ndims; ++i){
		if (i == 2){
			size_output.dims[i] = adjoint ? (ych = kch1) : (ych = kch2);
		}
		else{
			size_output.dims[i] = shm_input->dims[i];
		}
	}
	size_output.type = shm_input->type;

	shm_array_t *shm_output = 0;
	if (make_shm_array (output, &size_output, &shm_output)){
		perror(output);
		return -1;
	}

	cplxf *x = (cplxf*) shm_input->data.c;
	cplxf *k = (cplxf*) shm_kernel->data.c;
	cplxf *y = (cplxf*) shm_output->data.c;

	for (int b = 0; b < batch; ++b){
		spirit (y, adjoint, k, ks1, ks2, kch1, kch2, x, xs1, xs2, xch);

		y += xs1*xs2*ych;
		x += xs1*xs2*xch;
		k += ks1*ks2*kch1*kch2;
	}

	return 0;
}
